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Abstract 

Snow cover area affects snowmelt, soil moisture, evapotranspiration, and ultimately streamflow. 
For the Distributed Model Intercomparison Project - Phase 2 Western basins, we assimilate 
satellite-based fractional snow cover area (PSCA) from the Moderate Resolution Imaging 
Spectroradiometer, or MODIS, into the National Weather Service (NWS) SNOW-17 model. 
This model is coupled with the NWS Sacramento Heat Transfer (SAC-HT) model inside the 
National Aeronautics and Space Administration's (NASA) Land Information System. SNOW- 17 
computes PSCA from snow water equivalent (SWE) values using an areal depletion curve. Using 
a direct insertion, we assimilate fSCAs in two fully distributed ways: 1) we update the curve by 
attempting SWE preservation, and 2) we reconstruct SWEs using the curve. The preceding are 
refinements of an existing simple, conceptually-guided NWS algorithm. Satellite fSCA over 
dense forests inadequately accounts for below-canopy snow, degrading simulated streamflow 
upon assimilation during snowmelt. Accordingly, we implement a below-canopy allowance 
during assimilation. This simplistic allowance and direct insertion are found to be inadequate for 
improving calibrated results, still degrading them as mentioned above. However, for streamflow 
volume for the uncalibrated runs, we obtain: (1) substantial to major improvements (64-81 %) as 
a percentage of the control run residuals (or distance from observations), and (2) minor 
improvements (16-22 %) as a percentage of observed values. We highlight the need for detailed 
representations of canopy-snow optical radiative transfer processes in mountainous, dense forest 
regions if assimilation-based improvements are to be seen in calibrated runs over these areas. 
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1. Introduction 

Snow stores and releases water at different time scales, providing a significant portion of 
runoff and the human water supply, both globally (e.g., one-sixth of the population: Bartnett et 
al. [2005]) and particularly over snowmelt-dominated regions (e.g., half in the Western US: 
Serreze et al. [2001]). Snowpack-caused runoff delay has major relevance for hydrologic models 
simulating downstream surface and ground water availability for coupling with water use models 
(e.g., Schmid and Hanson [2009] coupling with irrigated agriculture model for Western US 
Sierra Nevada). 

In many current hydrologic and land surface models, mountain region simulations reflect 
snow parameterization limitations related to inadequate model representation of forcing, 
topography, net radiation, fractional snow cover and interactions with mainly forest vegetation, 
among other factors [e.g., Nijssen et al., 2003]. However, accurate spatial distribution of snow 
properties in these regions, including fractional snow covered area or extent (fSCA) is 
considered important for adequately predicting snowmelt, soil moisture, evapotranspiration, and 
ultimately streamflow. Besides the commonly studied snow depth, fSCA also directly influences 
snow energy and mass balances to provide melt and streamflow [Bios chi, 1991]. The fSCA also 
strongly affects the land surface energy balance in weather and climate models through high 
snow albedo. 

This study focuses on improving modeled streamflow by assimilating satellite-based 
fSCA for the challenging case of mountainous, dense forests. Past studies in such basins have 
struggled to demonstrate improvements in late melt-season streamflow simulations through late 
melt-season fSCA assimilation. Example workarounds accordingly employed are switching off 
the assimilation during the late melt season to avoid degrading simulations, and using 
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assimilation-based forecasts starting earlier in the melt season. With this in mind, we attempt to 
reject the null hypothesis that the state assimilation of fSCA during the late melt season does not 
improve late melt-season streamflow simulations in a prediction sense. Section 2 provides a 
literature review of the effect of densely forested canopy on remotely sensed snow extent. This 
is followed by a review of the influence of fSCA assimilation on modeled streamflow in Section 
3. Next, sections 4 and 5 detail the study basins and data, and the model and assimilation 
methodology used. Finally, sections 6 and 7 cover the analysis of results, discussions and 
recommendations for future work. 

2. Snow extent data and the densely forested canopy effect 

Satellite data have unmatched spatial and temporal coverage, providing fSCA estimates 
which could potentially improve hydrologic predictions. However, such fSCA maps can feature 
significant topography-caused noise over mountains [ Nagler et al., 2008]. Dense canopy 
coverage can have an even greater effect: existing simple representations of canopy radiative 
transfer in most remotely sensed snow product algorithms are inadequate for use over forests, 
wherein "remain one of the largest sources of uncertainty in the remote sensing of snow" for both 
electro-optical and microwave sensors [ Essery et al., 2009; Also Foster et al., 1991; Flail et al., 
1998; Simic et al., 2004; Vikhamar and Solberg, 2002 etc.]. Also, especially underneath the 
canopy, blockage from trees makes forest snow mapping more difficult at large off-nadir sensor 
view angles (e.g., for the Moderate Resolution Imaging Spectroradiometer or MODIS sensor; 
http://modis.gsfc.nasa.gov/) [e.g., Hall et al., 1998, 2001; Liu et al., 2008; Dozier et al, 2008]. 

Remotely sensed snow mapping uncertainty depends mainly on forest canopy type and 
density, e.g., microwave estimates degrade for densities above 60-70% [ Cline et al., 2004; 
Pulliainen et al., 2001], and electro-optical instruments like MODIS underestimate snow for 
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closed-canopy evergreen forests when compared against open-canopy areas [Simic et al., 2004]. 
Limited infonnation gleaned from viewable canopy gaps is insufficient to establish the below- 
canopy snow extent [e.g., Essery et al.; 2009]. Some factors causing differences between snow 
extents on and below the canopy include: canopies shading underlying snow from both direct 
and diffuse solar radiation [e.g., Pomeroy and Dion, 1996], and the nature of net radiation in 
dense canopies where sometimes increased longwave radiation [Pomeroy and Granger, 1997] 
and even convective terms [Sicart et al., 2004] potentially outweigh the effects of decreased 
shortwave radiation. 

Field studies have investigated the relative snow amounts, extents and durations present 
in open areas, forest canopy gaps, above the canopy, and below the canopy [see citations in the 
modeling study Rutter et al., 2009; Varhola at al., 2010]. While adjusting the remotely sensed 
snow depth using field measurements of forest cover has been attempted, e.g., for passive 
microwave observations [Armstrong and Brodzik, 2001; Clifford, 2010 etc.], the corresponding 
fSCA adjustment has received limited attention other than the simple assumption that the fSCA 
value in the forested portion of the pixel is same as that in the entire pixel [Molotch and 
Margulis, 2008]. Studies attempting better field-based fSCA adjustment include those by Liu et 
al. [2008] investigating the dependence of the sensor-viewable snow extent on the view angle, 
and Rice et al. [2010] conducting intensive in-situ gridded measurements of snow presence to 
adjust biases in the MODIS-derived MODSCAG product [Painter et al., 2009] for hydrologic 
modeling over the Sierra Nevada. 

Note that snow-vegetation radiative interactions can also complicate fSCA mapping 
indirectly based on snow-free vegetation properties and extents. For example, Robin et al. [2007] 
found vegetation index retrievals complicated by the presence of coniferous forest snow. 
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Efforts to improve snow parameterizations have traditionally been associated with open 
areas, and hence are mostly inadequate for forested areas [ Essety et al., 2009] where the above- 
mentioned knowledge gap relating overall snow extents to those in viewable canopy gaps 
remains. One option to obtain or properly utilize the remotely sensed fSCA is to use detailed 3- 
dimensional canopy radiative transfer \Pinty et al., 2004] techniques, such as a geometric 
optical-radiative transfer component using view angles, canopy profile and density, etc. to 
calculate variables like the viewable gap fraction [Liu et al., 2008], This can be applied either in 
the remotely sensed product algorithm or the corresponding post-processor, or in the hydrology 
model itself as in Hardy et al. [1997]. 

3. Existing studies of the influence of fSCA assimilation on simulated streamflow 

Most existing fSCA assimilation studies [e.g., Liston et al., 1999; Rodell and Houser, 
2004; Andreadis and Lettenmaier, 2006; Molotch and Margulis, 2008; Su et al., 2008; De 
Lannoy et al., 2012] focused on improving modeled snow water equivalent (SWE), with a 
minority focusing on improving streamflow. Non-assimilation studies focused on improving 
streamflow include those where ISC A was used directly as an input (forcing) [e.g., Li and 
Williams, 2008], or for model calibration [e.g., Udnaes et al, 2007; Parajka and Bldschl, 2008; 
Finger et al., 2011]. Below, we only discuss studies of the effect of fSCA assimilation on the 
modeled runoff/streamflow. 

Studies assimilating ISC A to improve streamflow had varying results: (1) Johansson et 
al. [2003] found ambiguous results when the model states are reset to those of the closest date 
where modeled fSCA equals the observation; (2) Huttunen et al. [2005] obtained mixed results 
on a boreal drainage basin; (3) Dressier et al. [2006] obtained reduced and increased streamflow 
simulation perfonnance respectively for relatively complex and relatively homogeneous terrains 
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in snow-dominated basins, with the performance of the fonner catching up with that of its non- 
assimilation counterpart when the assimilation was limited to the accumulation season; (4) Clark 
et al. [2006] found only minor improvements near the end of snowmelt season; (5) McGuire et 
al. [2006] showed the best improvements in forecasts made earliest in the melt season; (6) 
Udnaes et al. [2007] found that only some updates improved runoff, and most frequently at large 
SCA values; (7) Zaitchik and Rodell [2009] compared runoff as an 'unconditioned' sensitivity 
exercise (i.e., no comparison against actual observations) to show an increase and improvement 
in simulated snow volume, though their Figure 4 shows difficulty in simulating late-season SWE 
over the West Coast; (8) Roy et al. [2010] observed varying improvement in the streamflow 
Nash-Sutcliffe coefficient and the Root Mean Square Error (RMSE) using MODIS and NOAA 
IMS (National Oceanic and Atmospheric Administration Ice Mapping System) snow area data; 
(9) Tang and Lettenmaier [2010] found streamflow simulation errors were not necessarily 
reduced, with 2-week lead-time forecasts during snow ablation improving but not seasonal 
forecasts; (10) Third et al. [2011] found that discharge improved using the particle fdter 
assimilation method but deteriorated with the ensemble Kalman filter; and (11) Arsenault [2011] 
also performed an unconditioned sensitivity exercise at SNOTEL-coincident locations and 
demonstrated a reduction in snow mass and hence runoff during melt. Below, we highlight the 
similarities and differences between these studies and this paper. 

Our study basins’ physiography includes the double complexity of being mountainous 
and densely forested, with earlier studies mostly featuring one of these factors. The only 
successful study incorporating both seems to be by Roy et al. [2010] conducted over mostly 
forested basins (the Du Nord watershed is 96.2% forest-covered with 47.7% evergreen, and the 
Aux Ecorces basin is 83.5% forested with 12.3% evergreen). Somewhat puzzling is their 
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improvement in streamflow even though their Fig. 3 shows MODIS SCA decreasing relatively 
early to 0. Their methodology inserts a large 4 cm snow value into a snow-free model grid box 
when MODIS indicates snow but the model does not (see Section 5b: our study insertion is 
~5mm). The authors explain that this value possibly compensates for factors like the 
underestimation of winter precipitation; however, streamflow in their control simulation (which 
lacks such a compensation) does not follow the early decrease of the MODIS SCA. We speculate 
that a significant reason for their streamflow improvement is their large 4 cm snow insertion, 
since, per Simic et al. [2004]: (1) their Figure 8 shows Canadian evergreen and deciduous forests 
having a comparatively higher percentage of snow commission error than of omission, and (2) 
the 1 cm threshold is most representative of areal snow cover within the pixel. 

Regarding the transition to bare ground, Clark et al. [2006], Roy et al. [2010] and Udnaes 
et al. [2007] (some basins) with differing streamflow improvement results faced the problems of: 
(1) a substantial proportion of spring streamflow occurring before any bare ground is exposed, 
and (2) the transition from 100% snow-covered to snow-free conditions occurring fairly quickly. 
The rapid snow melt is evident in the Zaitchik and Rode/ 1 [2009] and Thirel et al [2011] studies. 

While other studies used fractional-valued spatial aggregates of binary snow/no-snow 
data (e.g., Roy et al. [2010], Thirel et al [2011]), our study assimilates fSCA at its base 
resolution. This is conceptually similar to Dressier et al. [2006], Udnaes et al. [2007] and 
Arsenault [2011], and to Clark et al. [2006] who used synthetic fSCA. Also, our snow model 
areal depletion curve (ADC) considers more depletion behavior characteristics than the simple 
ADC typically used, incorporating: (1) interannual consistency in the spatial pattern of relative 
snow amounts through a simple uniform scaling of the ADC shape across years, and (2) a 
temporary deviation of fSCA off the curve along a 'new snow' line to 100% and back when a 
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minimum new snowfall occurs on a partially bare area. While the ’new snow' line concept rarely 
exists in other studies (e.g., Huttunen et al. [2005]), the ADC shape preservation and scaling 
between years is present in the Johansson et al. [2003], Dressier et al. [2006], McGuire et al. 
[2006], and Tang and Lettenmaier [2010] studies (see Andreadis and Lettenmaier [2006] for 
details on the latter two). Typically, other models, e.g., the Noah LSM ( Chen et al. [1996]; Ek et 
al. [2003]) investigated by Zaitchik and Rodell [2009], have fixed ADCs that depend only on 
land cover (i.e., the depletion curve SWE at the maximum fSCA is a fixed maximum and does 
not reflect the variation in these maximum SWEs among different accumulation periods). 

Previous studies have also used multiple fSCA assimilation methods including: (1) 
resetting model states to those of the closest date where the modeled fSCA equals the 
observation (. Johansson et al. [2003]), (2) sequentially correcting the simulation by changing 
inputs such as temperature, precipitation and potential evaporation so that the observed SCA and 
other measured states agree with simulations (. Huttunen et al. [2005]), and (3) multiplying a 
snow depth with the snow extent to get effective SWE for a pixel as in Dressier et al. [2006]. 
Regardless, the assimilation approach is typically either threshold-based or ADC-based using the 
fSCA vs. SWE-based ADC. While our (and most other) studies follow the latter approach, some 
studies utilized rule-based updating where the snow presence decision is based on an fSCA 
threshold at the model’s resolution, e.g., McGuire et al. [2006], Roy et al. [2010] and Tang and 
Lettenmaier [2010]. 

Our study uses quality-controlled forcing data from Phase 2 of the Distributed Model 
Intercomparison Project (DMIP2; http://www.nws.noaa.gOv/oh/hrl/dmip/2/ ; [ Smith et al., 
2010a]), and, therefore, does not attempt to nudge meteorological forcings a day or more in 
advance of assimilation like Huttunen et al. [2005], Udnaes et al. [2007] and Zaitchik and Rodell 
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[2009]. We apply a below-canopy allowance to observed fSCA during assimilation to account 
for below-canopy snow, while other techniques include the pixel canopy correction factor for 
scaling up the snow extent by Dressier et al. [2006]. Like our study, most previous assimilation 
studies were not ensemble-based. The exceptions were Clark et al. [2006] and Thirel et al 
[2011], and also Arsenault [2011] which used an ensemble scheme in addition to direct insertion. 
Our study is a first-cut direct insertion attempt at using and updating the depletion curve and 
model states, with development of ensemble-based techniques left for future work. 

4. Basins and data used 

As mentioned above, our study is conducted over the DMIP2 Western US basins and is 
forced with hourly DMIP2 temperature and precipitation on the ~4 km Hydrologic Rainfall 
Analysis Project (HRAP: Greene and Hudlow [1982]) grid. We assimilate MODIS satellite- 
based fSCAs into this distributed DMIP2-type model simulation. In addition to comparisons of 
modeled versus observed streamflow, modeled SWE values are also compared against SNODAS 
SWE analyses [ Carroll et al, 2001; NOHRSC, 2004]. 

a. DMIP2 Western basins (North Fork American and the nested East Fork 
Carson) 

The DMIP2 hydrological model intercomparison experiment provided a framework to 
test many distributed models with operational quality data, for meeting the needs of operational 
National Weather Service (NWS) forecasters. Figure 1 shows the two hydrologically complex 
DMIP2 Western basins located in the Sierra Nevada Mountains. Their hydrologic complexities 
include snow, orographic precipitation, rain-snow partitioning, forest canopies, steep slopes and 
other complex terrain features. 
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Henceforth, the terms ’Carson’ and ’American’ refer respectively to the East Fork of the 
Carson River and the North Fork of the American River. The basins lie on either side of the 
Sierra divide [ Simpson et al., 2004], with the Carson on the eastern, leeward, rain-shadowed side 
draining south to north, and the American on the wetter, western, windward side draining 
westward. Average basin elevations are 2417 m and 1270 m respectively for the Carson and 
American, with respective ridgelines near 3400 m and 2700 m. Although situated close 
geographically, their hydrologic regimes differ greatly: the high-altitude Carson is snow- 
dominated, while the lower-elevation American has rain and mixed rain-snow events. Annual 
precipitation ranges from 560-1244 mm over Carson and from 813-1651 mm over the American. 

Both basins are geologically dominated by the granodiorites of central Sierra Nevada, 
with additions from volcanic and metasedimentary rocks for Carson and American respectively 
[Jeton et al., 1996]. The Carson has shallow sandy and clay soils, while American is 
characterized by clay loams and coarse sandy loams. Vegetation varies from sub-alpine and 
alpine conifer forests and meadows at the upper elevations of both basins, to chaparral-sagebrush 
rangelands in the lower elevations of the Carson, and pine-oak woodlands and shrub rangelands 
in the lower regions of the American. Based on the UMD vegetation continuous fields data 
( http://glcf.umiacs.umd.edu/data/vcf/) , the Carson’s higher elevation grid cells that receive most 
of the snow have forest canopy coverage ranging from 65 % to 70 %, and the lower elevations 
have coverage ranging from 30 % to 60 %. The American’s forest canopy coverage has a more 
uniform value of approximately 73 %. Both basins are largely unregulated [Jeton et al., 1996; 
Dettinger et al., 2004], although there do exist a few small reservoirs and diversions. 

While springtime snowmelt runoff dominates Carson streamflow overall, the largest 
peaks stem from rain-on-snow events. In contrast, about two-thirds of streamflow in the 
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American stems from wintertime rainfall and snowmelt runoff, and less than one-third from 
springtime snowmelt runoff ( Dettinger et al., 2004). Hourly streamflow in the Carson is 
measured by two nested USGS instantaneous stream stage gauges: the downstream Gardnerville 
gauge 10-30900 in Nevada and the upstream Markleeville gauge 10-308200 in California. These 
gauges drain, respectively, 922 and 715 km , while the American gauge 11-427000 at North 
Fork Dam 1 1-427000 in California drains 886 km 2 . 

b. The MODIS MOD10A1 fractional snow extent product 

The fractional snow cover area (fSCA) is drawn from Collection 5 MOD10A1 500 m 
daily fSCA [ http://modis-snow-ice.gsfc.nasa.gov ; Riggs et al, 2006] obtained by the MODIS 
electro-optical sensor on the Terra spacecraft. The MOD10A1 fSCA calculation involves a 
statistical-linear equation ( Salomonson and Appel [2004; 2006]) to the Normalized Difference 
Snow Index (NDSI). The assumed 'truth' data for this equation were binary snow presence values 
from the Landsat-7 Enhanced Thematic Mapper-Plus (ETM+; http://landsat.gsfc.nasa.gov/) . Like 
MODIS, ETM+ also does not properly sense snow under forest canopy for the Sierra Nevada, 
but its finer 30 m resolution can lead to slightly more accuracy than the MOD10A1 fSCA. 

Aggregation of the MOD10A1 data to a spatially coarser resolution (e.g. our ~4 km 
HRAP grid) can decrease confidence in the coarser estimates if any 500 m cloud covered pixels 
exist. In the 0.05° MOD10C1 MODIS fSCA product, the confidence index (Cl) and the snow 
extent are, respectively, the percentage of cloud-free pixels and binary snow pixels among all 
land pixels. MODlOAl's mutually exclusive mapping of cloud and snow pixels means the 
MOD10C1 extent calculation assumes no snow under cloud pixels, whereas a possibly better 
assumption is that the proportion of snow-free land versus snow-covered land among the cloud- 
free pixels is also maintained among the cloud pixels. We accordingly calculate the Cl but by 
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using (EC A instead of binary snow values, where the coarser grid’s fSCA is the average of the 
cloud-free pixels' (EGAs. The Land Information System software (LIS: Kumar et al. [2006]; 
Peters-Lidard et al. [2007]) directly reads the MOD10A1 via our additional coded internal 
reader utilizing HDF-EOS (Hierarchical Data Format - Earth Observing System) library tools. 
The LIS integrates satellite- and ground-based observational data with advanced land modeling 
techniques to produce optimal fields of land surface states and fluxes. Our study assimilates 
fSCA values with a high Cl value of 80% (compared to values of 50-80% in Andreadis and 
Lettenmaier [2006], 90% in Hall et al. [2010]). Lower Cl thresholds have been used in other 
studies in an attempt to retain more information during snowfall when clouds are present (e.g., 
6% by Rodell and Houser [2004] and Zaitchik and Rodell [2009]). 

c. The SNODAS snow analysis product 

The SNODAS [ Carroll et al., 2001; NOHRSC, 2004] is a spatially distributed modeling 
and data assimilation framework which attempts to provide the best possible physically 
consistent (i.e., snow energy- and mass-balance based) estimates of snow properties on a 1 - km 
grid throughout the US and parts of Canada to support hydrologic modeling and analysis. It 
assimilates observations from ground stations, satellites, and airborne passive gamma-ray 
sensors. These gamma-ray sensors may potentially sense snow better under the canopy in 
mountainous forests than do their electro-optical remote sensing counterparts (e.g., MODIS). 
Model forcings includes downscaled fields from numerical weather prediction models, surface 
weather observations, satellite-derived solar radiation and radar-derived precipitation. 

5. Model and Assimilation Methodology 

a. NOAA OHD’s hydrology component models and their porting into LIS 
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The National Oceanic and Atmospheric Administration’s Hydrology Laboratory Research 
Distributed Hydrologic Model (HL-RDHM) is an updated version of the HL-RMS (Hydrology 
Laboratory-Research Modeling System: Koren et al. [2004]). It includes gridded hydrological 
model features and components. To exploit LIS data handling and assimilation capabilities, we 
ported into LIS the following HL-RDHM components: the gridded forcing reader, the SNOW- 17 
snow model [Anderson, 1973, 2006], and the Sacramento Heat Transfer (SAC-HT: Koren et al. 
[2007]; Burnash [1994]) surface water balance model. The hill- and channel-slope runoff routing 
post-processor utilizes a physically-based kinematic wave approach that provides unconditional 
stability [Reed, 2003; Koren et al., 2004]. 

Meteorological inputs are first passed to the SNOW- 17 model which represents the 
snowpack as a relatively simple one-dimensional bulk layer with an associated water holding 
capacity. Rain-on-snow periods use energy balance-based melt computations, while non-rain (or 
negligible rain) periods utilize a temperature index equation. During sub-freezing temperatures, a 
heat deficit function tracks the energy needed to return to an isothermal, or ripe, snowpack state. 
This deficit is a function of the temperature gradient between the snow surface temperature (the 
minimum of the air temperature and 0 °C) and a computed snowpack antecedent temperature 
index. Liquid water storage calculations using the heat deficit ensure preferential refreezing of 
any surface melt or rain before adding to the liquid water storage and outflow (in contrast to 
Slater and Clark [2006] where contribution to liquid water storage or outflow is preferred to 
refreezing within the snow). Any excess water is lagged in time and attenuated during 
transmission through the snowpack according to the ratio of the ice water equivalent and this 
excess water. SNOW- 17 has 22 parameters (Table 1). The shape of the areal depletion curve 
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(ADC) is described by 1 1 of these parameters, and the ADC and its adjustment during fSCA 
assimilation are described in a following subsection. 

Snowmelt from SNOW- 17 and precipitation over snow-free areas force the SAC-HT 
water balance model. SAC-HT has upper and lower zone tension and free water storages 
representing a relatively thin upper layer, and a much thicker lower layer. The model simulates 
evaporation, saturation-excess infiltration, percolation, fast runoff response components (surface 
runoff and direct runoff from impervious surfaces), slower runoff responses (interflow, 
supplemental and primary baseflow), and non-channel subsurface outflow. 

A grid cell's drainage density parameter delineates conceptual overland flow plane 
hillslopes. SAC-HT hillslope fast response runoff drains into a conceptual channel having the 
cell diagonal as its maximum length. Hillslope slow response runoff enters the channel directly, 
b. The SNOW-17 areal depletion curve and direct insertion assimilation of fSCA 
Different areal depletion curve (ADC) types relate the fSCA to the current fraction of 
total seasonal runoff, the time of the year, degree-days, and the SWE. Common operational 
ADCs as used in SNOW-17 calculate fSCA from the modeled SWE (i.e., fSCA is diagnostic ). 

For a given area, the general shape of such fSCA vs. SWE ADCs is generally similar 
across years, reflecting the fairly high interannual consistency in the spatial variations of the 
relative snow amounts [e.g., Liston, 1999, 2004; Luce et al, 1999; Luce et al., 2004; Anderson, 
2006; Kolberg and Gottschalk, 2010; Homan et al, 2011]. Many factors affecting accumulation 
and melt combine towards this year-to-year consistency, including temperature, storm type and 
direction, wind speed, cloudiness, dew-point temperature, topographical factors like elevation, 
slope and aspect, and vegetation [e.g., Konig and Sturm, 1998; Luce et al., 1998; Deems et al., 
2008; Griinewald et al, 2010; Rice et al., 2010; detailed review in Clark et al., 2011]. In the 


Page 16 of 54 


SNOW- 17 ADC, the Snow Index (SI) parameter is the upper limit of the Areal Index A;, which 
is itself the variable minimum seasonal SWE accumulation at 100% fSCA. So SI is the minimum 
areally-averaged SWE above which fSCA is always at its maximum possible value (100 %), and 
below which the fSCA is 100 % or less depending on SWE and Aj [Anderson, 2006] (see SI limit 
in Figure 2a). At maximum fSCA, when Aj updates upward during the control run (i.e., non- 
assimilation) accumulation, the entire ADC instantaneously steepens. At non-maximum fSCA, 
the ADC has an additional ’new snow' trajectory behavior when fSCA temporarily becomes 
100% from an arbitrarily set 0.2 mm/hr minimum snowfall on a partially bare area. This 100% 
fSCA remains until an arbitrarily chosen value of 25% of this newly fallen snow melts, and then 
the trajectory linearly returns to the earlier pre-new snow trajectory point on the ADC (Figure 
2a). 

Our direct insertion technique is based on the NWS operational manual modification 
(MOD) of state variables [Anderson, 2002; 2006] called the Areal Extent of Snow Cover Change 
(.AESCCHNG). Note that this MOD notation differs from the MODIS MOD10A1 data. The 
.AESCCHNG preserves SWE by assuming no intrinsic SWE information in the observed fSCA. 
Assimilating fSCA gives different ADC trajectory alterations cases depending on observed fSCA 
relative to modeled values. These alterations change snow processes, SWE and timing of melt 
down the line, eventually impacting streamflow. The case depicted in Figure 2b shows the pre- 
assimilation (i.e., “Old”) SNOW-17 fSCA where the entire ADC is instantaneously able to 
steepen and flatten during both snow accumulation and melt phases. By comparison, the control 
run ADC continuously steepens and only during accumulation. Thus the observed fSCA instead 
of the A; controls the trajectory change, i.e, fSCA is now prognostic as compared to being 
diagnostic in the control run. Figures 2c-d depict cases starting on the 'new snow' line where 
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.AESCCHNG adjusts either this line, or the ADC (i.e. new A;), depending on whether the 
assimilated fSCA is respectively greater or less than the fSCA corresponding to the SWE on the 
pre-assimilation ADC (i.e., see the dashed ’Threshold fSCA’ line in Figure 2 legend). 

Note that .AESCCHNG is used only in cases where both the modeled and observed fSCA 
are non-zero (i.e., mutually exclusive snow/no-snow cases between model and observations). 
Additionally, .AESCCHNG could allow a new A; beyond the realistic upper limit of SI. In 
response to this, we simultaneously enhance .AESCCHNG to avoid this Aj upper limit problem 
and the inaction during the mutually exclusive snow/no-snow cases. This is accomplished by 
implementing an SI constraint when the SWE goes above SI to enforce the post-assimilation 
fSCA at 100%, irrespective of the observed ISC A (Figures 2e-f). Similar to Rodell and Houser 
[2004] and Zaitchik and Rodell [2009], we handle the mutually exclusive cases of snow/no-snow 
through allowing ’unavoidable’ alterations to the SWE (not shown): [1] For modeled extents at 
zero, we add SWE derived from non-zero observed (SC A and an ADC defined by A, at 5 mm, 
and [2] For observed fSCA at zero, we remove existing modeled SWE. 

In addition to these .AESCCHNG-based runs that assume no intrinsic SWE information 
in the observed fSCA, we also test the opposite assumption that combines fSCA with the existing 
modeled ADC to reconstruct SWE (Figures 2g-h). This differs from the operational Water 
Equivalent Change (.WECHNG) MOD that adjusts the ADC using observed SWE. 

c. Hydrological model parameterization, including the areal depletion curve 
Tables 1-3 lists the parameter names, corresponding values, ranges and scalar multipliers 
for the uncalibrated ( a priori ) and calibrated runs. The Carson, with two nested streamflow 
gauges (Figure 1), serves as the domain for two uncalibrated and three calibrated control runs. 
The extra calibrated control run makes use of the interior Markleeville gauge as a ’blind’ 
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streamflow simulation point testing the Gardnerville domain calibrated parameters. The single 
American gauge has only one uncalibrated and one calibrated control run. The hitherto 
undetailed process of creating the uncalibrated and calibrated parameters occurred as part of the 
DMIP2 Western basins study ( Smith et al, 2010b), so is briefly outlined below for completeness. 

The SNOW- 17 a priori values for MFMAX, MFMIN and MBASE are obtained by 
combining a derived energy-based temperature index snowmelt equation for spring snow pack 
with the SNOW- 17 temperature index model, using available spatial physiographic data, and 
wind speed climate grids from the North American Regional Reanalysis (NARR) [ Mizukami and 
Koren, 2008]. With the exception of the ADC, the remaining parameters were derived from 
lumped-calibrated values from the California Nevada River Forecast Center and knowledge of 
the local snowfall and wind climatologies. The SAC-HT a priori physically-based parameters 
are derived from physiographic properties like soil texture and land cover [ Koren et al, 2000; 
Zhang et al, 2011]. The routing a priori parameters are from DEM, land use and empirical 
equations based on channel hydraulic data [Koren et al, 2000], 

In the manual American calibration process, the distributed a priori values are first scaled 
to match pre-existing calibrated lumped values. All scalar parameters are then adjusted to 
replicate observed hourly stream flow using the following sequential methodology for lumped 
models [ Smith et al, 2003]: Removal of major biases and errors, matching of SAC-HT base 
flows, general SNOW- 17 calibration, further calibration to match in-situ observed SWE, 
matching of the faster responding surface runoff, and adjustment of routing parameters. 
Streamflow matching is done using visual comparison and multiple goodness-of-fit measures at 
different stages, e.g., base flow statistics, seasonal statistics, flow interval biases, run period 
accumulated error, overall bias and the modified correlation coefficient [McCuen and Snyder, 
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1975]. The Carson calibration is similar, except: (1) Parameters from two elevation zones in a 
preexisting lumped calibration give two corresponding zone fields (though single multipliers 
then varied both zones together), and (2) A mixture of manual and Simplified Line Search 
automated approaches [Kuzmin et a/., 2008] provides the final calibrated distributed values. 

A principal difference between our study's and DMIP2's uncalibrated Carson 
configuration is our non-zero SI value (versus the DMIP2 a priori SI value of zero). This change 
was essential, as a zero SI means only 0% or 100% values of fSCA are allowed, effectively 
reducing the fSCA information during assimilation to rule-based binary values as in Rodell and 
Houser [2004]. We instead use the DMIP2 calibrated SI (and associated ADC) for both 
calibrated and uncalibrated Carson runs. The DMIP2 uncalibrated zero SI is simplistic and 
reflects typical operational constraints on the assignment of SI and ADC. More broadly, where 
skill level or time available is not optimal, a simplistic SI assignment can be either (1) a 0 as 
mentioned above where the basin-wide fSCA becomes the percent of pixels with snow presence, 
or (2) a very high value like 9999 or 999 mm (e.g. 999 for American) that is always above 
existing SWE values so that bare ground is exposed immediately as melt starts. This second 
option is useful when future plans include refined SI calibration using multiple years of high- 
snow data. Snowmelt and streamflow can be very sensitive to SI (. Anderson [2002]). For 
example, a low (high) SI gives corresponding lower (higher) Aj, higher (lower) fSCA derived off 
the depletion curve, higher (lower) snowmelt and higher (lower) streamflow. 

Similar to the simplistic SI assignment mentioned above, all American runs feature a 
simplistic straight line ADC assignment [Anderson, 2002], which perfonned well in DMIP2. 
However, for the Carson, the more sophisticated DMIP2-calibrated ADC was derived using 
previously established ADCs for well-studied, calibrated basins in the same region, in addition to 
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information on terrain, vegetation, climatic conditions, satellite measurements, and any available 
ISC A (e.g., from the National Operational Hydrologic Remote Sensing Center or NOHRSC). 

d. Snow assimilation allowance under dense forest canopy during snow melt 

The top panel of Figure 3 shows a bird’s-eye nadir view of satellite-sensed fSCA 
(fSCA S ateiiite) components in a conceptual grid square area, based on vertically projecting forest 
canopy and snow extents onto the ground surface: areas I, II, IV and III denote no-canopy no- 
snow, no-canopy snow, canopy no-snow, and canopy snow respectively. Values IV and III are 
based above the canopy. The total no-canopy and canopy areas are respectively (I+II) and 
(IV+III). Satellite sensors miss the below-canopy snow component V of the true fSCA (fSCA true ) 
in the lower panel of Figure 3. As such, fSCA sat eiiite is the sum of only II and III (II and III values 
are not separately given). During accumulation (Figure 3A), snow covers almost the entire grid 
square (i..e, I & IV are 0 or negligible), with the above-canopy amount (III) usually greater than 
the below-canopy amount (V), indirectly making fSCA sa temie equal to fSCA true . Assimilating 
fSCA in this stage is done without any melt time allowance, contrary to what is done during the 
snow melt phase described below. 

We assume that the late melt stage roughly occurs when the above-canopy snow III has 
become less than the below-canopy V (Figure 3C). Since the temporal variation of II, III and V 
values are unknown under existing model and data constraints, we provide an unavoidably wide 
melt-period upper allowance equal to the forest canopy area on the fSCA sate iiite to be assimilated 
(hereforth denoted the ’UnsensedSnow' effect). Note that the ideal allowance should be a reduced 
value of IV (or canopy area minus the III value) instead of the entire canopy area, but this 
reduction is unachievable since the III value is not separately known. The assimilation is limited 
to constraining post-assimilation fSCA to outside this allowance-included range bounded by 
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fSCA sa teiiite and (fSCA sa teiiite + canopy area). So when pre-assimilation modeled fSCA (fSCA mo d e i) 
is greater than the (fSCAsateiiite + canopy area) upper bound, the assimilation process reduces it to 
equality with that bound: the assimilation inaccuracy is the entire canopy area minus the 
difference between above-canopy III and below-canopy V. And when the pre-assimilation 
fSCA mo dei is less than the fSCA sate iiite lower bound, the assimilation process increases it to the 
fSCAsateiiite value: the assimilation inaccuracy is the difference between III and V. 

During early melt periods assumed to occur when the above-canopy snow III is greater 
than the below-canopy V (Figure 3B), an accumulation period-type assimilation without 
allowance seems suitable. However, since the timing of exact transition between early and late 
melt is not known under current model and data constraints (or when the relative dominance 
between III and V switches), we again implement a late melt period-style allowance. So when 
pre-assimilation fSCA mo d e i is greater than the (fSCA sa teiiite + canopy area) upper bound and the 
assimilation process reduces it to equality with that bound, the assimilation inaccuracy increases: 
it is now the canopy area plus the difference between above-canopy III and below-canopy V. For 
the case where the pre-assimilation f$CA mo dei is less than the fSCA sate iiite lower bound, direct 
insertion assimilation inaccuracy is absent. Note that this entire first-cut approach violates a 
SAC-HT evapotranspiration demand calculation assumption that the snow areal proportion is the 
same in the canopy and no-canopy areas, a shortcoming that future refinements to the 
assimilation technique will address. 

e. Description of the direct-insertion assimilation simulations 
All simulations span September 1995 to September 2006, mirroring the DMIP2 
simulation period. Data assimilation starts when MODIS fSCA becomes available (March 2000 
for Terra MODIS), and we stress that a long spinup period is not required for our application. 
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Thus six complete water years (WYs: Octobers to Septembers) are available for evaluating the 
assimilation procedure (WYs 2001-06). The SNODAS data series commences in October 2003, 
providing only three WYs for SWE comparison (WYs 2004-06). Assimilation results for the 
SNODAS 3-WY period subset are similar to the 6-WY values and not reported. Table 4 lists the 
multiple runs. The control run corresponds to Figure 2a, the Mod assimilation run A to Figures 
2b-d, the enhanced Mod runs B and C to Figures 2e-f, and the SWE-recreation assimilation run 
D to Figures 2g-h. The uncalibrated control run serves as the baseline against which to judge the 
performance of the uncalibrated assimilation runs' suite, and of an extra calibrated control run 
(denoted by Z) against which this suite can be compared. Note that runs B and C differ from A in 
having addressed the A; upper limit and snow/no-snow issues (see Section 5b). 

f. Assimilation performance evaluation measures 

We compute event-based and continuous measures. Events considered are an evaluation 
period subset of DMIP2 events. We adapt DMIP’s event-based improvement measures [Smith et 
al., 2004] to calculate streamflow error improvement from assimilation, against the control runs 
(runs here are spatially distributed). These measures include: flood runoff improvement, peak 
flow improvement and peak time improvement (Appendix A equations 1-3). We also construct 
similar additional performance measures where improvement is calculated as a percent of the 
control run residual values instead of the observation values (Appendix equations 4-6): we 
henceforth specifically call these 'residual improvement measures' as compared to the former 
ones called 'improvement measures'. These residual improvement measures are more relevant for 
determining how much the deviation of the control values from the observations is reduced by 
the assimilation, especially where these control value deviations, and so improvement measure 
values are low. In other words, where the residual improvement measure is a high y %, a low 
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improvement measure at x % can theoretically improve to a value of only around x*(100/y) % 
when the assimilation results exactly match observations. This x*(100/y) % may still be of the 
same order as x %, indicating little available room for improvement. 

Time-continuous streamflow improvement measures for a multiple Water Year (WY) 
period are computed as the average of corresponding single WY values. So streamflow bias 
improvement and streamflow bias residual improvement are the respective equivalents of the 
above-mentioned event-based measures concerning flood runoff (Appendix equations 1 and 4). 
We also calculate the assimilation run's modified correlation coefficient and its improvement 
over the control run (equations 7-8). Next, we calculate the percent of time steps where the 
assimilation run is better than the control run, and the streamflow residual improvement that 
occurs within those time steps (equations 9-10). 

We also use time-continuous SWE improvement measures to assess basin-mean SWE. 
These are the SWE bias improvement, SWE bias residual improvement, SWE centroid date 
improvement and SWE centroid date residual improvement (equations 11-14). Note that possible 
limitations of SNODAS SWE (very few gamma-ray flights per melt season etc.) can mean that 
conclusions based on the SWE measures' values are possibly 'weaker' and so used more for 
comparison than for evaluation. 

6. Analysis of results 

Figure 4 shows example melt season time series for the uncalibrated Carson. Run B SWE 
and streamflow decrease much faster (Figures 4b-c), due to MODIS IlSCA values decreasing to 
0. This behavior is irrespective of calibration or basin, and clearly shows the importance of 
allowance for below-canopy snow in satellite-based fSCA assimilation during melt. Figures 4d-f 



Page 24 of 54 


plot monthly mean streamflow residual improvement against observed streamflow, confirming 
that the crucial melt season has well-performing assimilation time steps for some runs. 

We do not list peak timing measure values: though residual improvements (RI t ) can be 
substantial, the corresponding actual improvements (I t ) are too low to be meaningful. This means 
that the observed timing is already well-reproduced here, even for the uncalibrated control. 
Values for all other measures are in Table 5 and Figure 5. Immediately noticeable is that while 
some uncalibrated run streamflows do show improvement, the calibrated run streamflows mostly 
degrade, sometimes by a large amount. The unavoidable allowance-related inaccuracy in our 
assimilation compromises the combination of detailed physics representations and accurate 
satellite fSCA (in addition to a well-calibrated ADC) that is required to improve streamflow. For 
the calibrated Carson example, this inaccuracy is likely larger than the difference between the 
pre-adjustment fSCA m<K i c i and the fSCA tru e, leading to degraded streamflow. The opposite is true 
for the well-performing uncalibrated Carson runs. We assume an improvement as significant or 
substantial if the percentage value for observation-related improvement measures (e.g., I y , not the 
residual improvement measure RI y in Appendix 1) is above five. This criterion cannot assign 
significance to the SWE centroid date measures since I c is in units of days. Hereafter, we discuss 
uncalibrated run results relevant to ungauged basins. 

The uncalibrated Carson run A (’MODIS: Mod') actually perfonns better than the control 
for SWE-related RIw, Iw, RIc and Ic as depicted in Table 5, but displays consistent degradation 
against the control run for all streamflow-related measures RI y , I y , Ir mo d, and r mo d, and low B t 
values. This apparent discrepancy is explained through visual examination of the high-biased 
SWE in uncalibrated run A which fails to completely melt off during the summer (not shown): 
since SWE in the control and in other assimilation runs is mostly under-biased against SNODAS, 
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higher SWE values from assimilation can be closer to observations and give a false impression 
of performance improvement (if streamflow scores are ignored). 

SWE improvement measures RIw, Iw, RIc and Ic in Table 5 consistently improve from 
run C ('MODIS: UnsensedSnow') to D ('MODIS To SWE: UnsensedSnow’) regardless of basin 
or parameter calibration (not ADC calibration). This indicates that fSCA sat eiiite recreates SWE 
instantaneously using the ADC better than through long-term ADC trajectory changes. 

The only assimilation run where a significant and consistent improvement occurs in both 
SWE and streamflow is the best-perfonning uncalibrated Carson run D (’MODIS To SWE: 
UnsensedSnow’). Corresponding SWE improvement measures RIw and Iw are 45 % and 9 % 
respectively (i.e., the average of the italicized ranges 42-48 % and 8-10 % in Table 5). Similarly, 
Table 5 run D values for RIc, Ic, streamflow RI y and I y are respectively 28 % (from 27-29 %), 
3.5 days (from 3. 4-3. 7 days), 68 % (from 64-72 %) and 17 % (from 16-18 %). Although the 
flood I y and I p values look significant for Gardnerville, their values for the SNODAS-available 3- 
year sub-period fall below the arbitrary 5 % threshold. The uncalibrated Carson run D also does 
not degrade performance in terms of correlation (Figure 5), and features a high percentage of 
time when its streamflow is better than the control run (B t ~70%). Run D likely performs better 
than runs A-C for the snow-dominated Carson because the below-canopy snow allowance 
supports better updates that instantaneously use the well-calibrated ADC to promptly propagate 
(iSCA values into better SWE values. This leads to better runoff and streamflow (dominated by 
springtime snowmelt runoff) from the large snowpack. Note that the performance of this 
uncalibrated Carson run D is close to that of the calibrated run Z across all improvement 


measures. 



Page 26 of 54 


Unlike the Carson wherein uncalibrated run D was the consistent best performer, 
attempts to select a similar run for the American are inconclusive. Uncalibrated run D is 
generally best in terms of SWE bias (Table 5 highlighted values are 25 % and 10 % for RIw and 
Iw respectively). However, uncalibrated run B is best in tenns of streamflow. Corresponding 
streamflow bias measures RI y and I y show strong improvement with respective values of 8 1 % 
and 22 %. By contrast, flood runoff and flood peak improvement values are modest: flood RI y , 
flood I y , Rip and I p are 8 %, 6 %, 6 % and 8 % respectively. The mixed rain-on-snow regime 
which characterizes the American leads to springtime snowmelt from snowpacks that are thinner 
than those of the Carson, and that account for less than one-third of the yearly runoff. Two-thirds 
of the annual runoff occurs earlier in the year from wintertime rainfall, where the assimilation 
allowance-related uncertainty is much higher than that in the Carson where the dominant runoff 
occurs towards the late melt season (refer to Section 5d). The difference in perfonnance ranks 
when judging against SWE and streamflow may indicate the presence of SNODAS SWE bias in 
such mixed rain-on-snow basins. 

7, Discussion and recommendations for future work 

Improving modeled streamflow through assimilation of satellite-measured snow area is 
an important research contribution [Andreadis and Lettenmaier, 2006; Simpson et a/., 2004]. We 
attempt this for the hydrologically challenging example of two Western US mountainous, 
densely forested basins (one snow-dominated and one rain- and snow-dominated) used in 
DMIP2. We apply a cautious below-canopy fSCA allowance during fSCA assimilation to avoid 
degrading simulated streamflow due to the lack of below-canopy measurements in the satellite 
fSCA product. The presence or absence of this allowance is one factor which differentiates the 
suite of direct-insertion assimilation runs presented in this study. 
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Compared to the impacts of SWE assimilation noted in other studies, the information 
content in observed fSCA can be generally modest towards improving streamflow simulations, 
especially over densely forested mountainous domains. Our study found fSCA assimilation 
mostly degraded streamflow for already well-performing calibrated runs where an accurate 
combination of model and data is required. This is because the below-canopy allowance 
inherently introduces an inaccuracy into the simulation that is likely more than the distance of 
pre-adjustment modeled fSCA (fiSCA mo d e i) from the true value (fSCA tme ). However, for 
uncalibrated runs, though the SNOW- 1 7/S AC-HT control run produces reasonably good 
streamflow perfonnance, fSCA assimilation still leads to: (1) substantial to major improvements 
(64-8 1 %) in streamflow volume as a percentage of the control run residuals (or distance from 
observations), and (2) minor improvements (16-22 %) in streamflow volume as a percentage of 
observed values. The apparent difficulty in achieving multi-objective (i.e. all aspects of 
streamflow and SWE) improvement is similar to that seen in earlier studies such as Third et al 
[2011] where application of an ensemble Kalman filter (EnKF) led to improved timing of snow 
extent depletion but degraded streamflow. These challenges notwithstanding, the direct insertion 
technique assessed in this study has the potential to benefit applications over ungauged basins 
(coming close to the impact that calibration has on model performance), especially for snow- 
dominated basins where significant streamflow occurs during the spring. 

Comparing this study's results against those of earlier studies is generally hampered by 
those studies having either non-no nnalized objective functions like root mean square error, or 
inappropriate nonnalized objective functions (e.g., Roy et al. [2010] use the Nash coefficient 
with the oft-used mean-of-observations benchmark rather than with more appropriate 
benchmarks like calendar day or climatology as explained by Sdiafeli and Gupta [2007]). Where 
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comparable, our results are at least of the same order. For example, Third et al [2011] employ 
both the ensemble Kalman filter (EnKF) and particle filter techniques on flatter basins, and their 
Table 3 MODIS assimilation results show the discharge score ratio bias changing from 0.13 for 
the control run to 0.09 for the assimilation. This means that their streamflow volume residual 
bias improvement (RI y ) is 100*(0.13-0.09)/0.13 = 31% and their streamflow volume bias 
improvement (I y ) is 100*(0. 13-0.09) = 4%. Our corresponding values for the uncalibrated Carson 
and American simulations are 68 % and 81 % respectively for RI y , and 17 % and 22 % 
respectively for I y (see respective italicized values of 72-64 %, 81 %, 16-18 % and 22 % in Table 
5). We also re-ran all the simulations reported in this paper with a lower (i.e., worse) confidence 
index acceptance threshold of 6 % for the MODIS data ( Rodell and Houser [2004]; Zaitchik and 
Rodell [2009]) and found that the conclusions remain valid with only a negligible degradation in 
the results. 

In practice, deciding on an appropriate run configuration based on our study suite will 
depend on the basin’s characteristics such as whether it features dense or sparse forest canopy, 
whether the majority of streamflow occurs during wintertime or springtime etc. Consideration of, 
and adjustment for, additional factors such as the time of the year, the ongoing stage of snow 
season accumulation or depletion, and whether the basin is snow-dominated or has mixed rain- 
on-snow events may also prove beneficial. For example, the American basin features significant 
wintertime streamflow, possibly coinciding with times when the intercepted snow extent and the 
amount on the canopy are greater than the below-canopy amount. In this situation, it could be 
advantageous to have the assimilation procedure switch off the below-canopy snow allowance 
during wintertime streamflow/snow depletion, but activate it during springtime snow depletion. 
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This allowance for prolonged below-canopy snow (e.g., Musselman et al. [2008]), missed by 
satellite instruments, is crucial when modeling the springtime melt period (Figure 5 b-c). 

Additionally, while the accurate determination of the timing of the movement of snow 
from the canopy downward to the ground would be vital in the selection of the proper 
assimilation method, our snow model (and most current hydrology models) do not include 
detailed canopy, snow and optical radiative transfer physics. We believe that such detailed 
process-oriented parameterizations providing infonnation on the snow extent above and below 
the canopy are the way forward for improving real-time streamflow simulations. It is possible 
that inclusion of such processes would negate the need for calibrating the areal depletion curve. 

In addition to more detailed process representation, additional ways to improve 
assimilation results involve the data (satellite-derived MODIS). The overall quality of ingested 
snow data can be increased by switching off the assimilation of lower satellite-observed values 
of fSCA (like those below 25% that typically have errors, Riggs et al. [2006]). Another way to 
use the distributed data connected to the assimilation technique is through joint spatial 
assimilation of the distributed extent values. Kolberg et al. [2006] reported on such 
improvements judged against fully distributed assimilation by Kolberg and Gottschalk [2006]. 
Joint spatial assimilation is possible and is a focus of ongoing work for the model used in our 
study. In such over-constrained problem involving the pixel and basin values, we are 
implementing a SWE equilibrium achievement between the individual and joint snow extents, 
SWEs and the depletion curve for robustness and hydrologic consistency. 

Assimilation results here could be further improved using better techniques to measure 
and differentiate precipitation into rain and snow, e.g., the vertically pointing radar-estimated 
bright-band height ( Lundquist et al. [2008]). Finally, better results are possible for hydrologically 
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less complex basins (the vast majority of basins, which feature gentler slopes or less forest 
canopy), given that improvements were possible even in our challenging scenario study. 


Appendix A 

1. [a] Flood runoff improvement and [b] Streamflow bias improvement (I y , %: streamflow 
improvement measure): 
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where Y 0J , Y C I and Y aji are the runoff volumes of observations, control run and assimilation 
run respectively for the z' th flood or Water Year (WY), Y <xavg is the average observed flood 
event runoff volume for N flood events or WYs. 

2. Peak flow improvement (I p , %: streamflow event-based measure): 
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where Q po ,i, Q pc ,i and Q paj i are the peak discharges of observations, control run and 
assimilation run respectively for the z' th event, Q p0j avg is the average observed peak discharge 
for N events. 

3. Peak time improvement (I t , hrs: streamflow event-based measure): 



where T po ,i , T pcj and T paj i are the time (hrs) of observations, control run and assimilation run 
respectively for the z th peak. 
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4. [a] Flood runoff residual improvement and [b] Streamflow bias residual improvement (RI y , 
% : streamflow measure): 
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5. Peak flow residual improvement (RI P , %: streamflow event-based measure): 
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6. Peak time residual improvement (RI t , %: streamflow event-based measure): 
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7. Modified Correlation Coefficient (r mod : streamflow continuous measure) 

From McCuen and Snyder [1975], this eliminates the regular correlation coefficient's 
tendency to be overly influenced by outliers and to be insensitive to differences in the 
hydrograph sizes: 
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where o siniJ and o 0 b s ,i are standard deviations of simulation and observation series 
respectively for the 7 th WY, p is corresponding regular correlation coefficient. 

8. Improvement in Modified Correlation Coefficient (Ir mo d: streamflow continuous measure) 
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where r mo d, c ,i and r mo d, a ,i are the control and assimilation run r mo d respectively for the z' th WY. 

9. Times where assimilation better than control (B t , %: streamflow continuous measure) 

.N b 


1 ^ N Bi 

B .= — Y(— ^-).100 


where Ns,i is the number of time steps that assimilation residual ( Q 0 /Y - Q a fi ) < control 
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residual ( Q o ji - Q c ji ) for the z' th WY, N, is corresponding number of time steps, Q 0 ji, Q 
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and Q a ji are the observed, control and assimilation run streamflows respectively at the j th time 
step. 

10. Streamflow residual improvement during times with better assimilation (RIq, %: streamflow 
continuous measure) 
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where Q 0 ji, Q c ,ji and Q a ,ji are the observed, control and assimilation run streamflows 
respectively at the j th time step of the z' th WY. 

1 1 . SWE bias improvement (I w , %: SWE continuous measure): 


I 


w 



N-W, 


om,avg 


•100 


where W omj , W cm ,i and W am ,i are the mean SWEs of reference (basin SNODAS), control run 
and assimilation run respectively for the z' th WY, W om , a vg is average reference mean SWE 
(basin SNODAS) for N WYs. 

12. SWE bias residual improvement (RIw, %: SWE continuous measure): 
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13. SWE centroid date improvement (Ic, days: SWE continuous measure): 
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where SCD 0 I , SCD cl and SCD aji are the SWE centroid date (days) of observations, control 
run and assimilation run respectively for the 7 th WY, 


SCD for the 7 th WY is 


Zb‘0 

i ,= i 

24.| 


ji is the j th hour, and W h , is the SWE at j,. 

14. SWE centroid date residual improvement (RIc, %: SWE continuous measure): 


RI C = 


Z |.SCD„,, - SCD CJ | - 1 SCD„, - SCD a i |) 

Z=1 
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1. Overview of the DMIP2 Western basins used in this study (adapted from Jet on et al., 1996 
and Smith et al, 2010b) showing stream gauge locations. 

2. Graphical depiction of direct- insertion update cases using example ADCs and ’new snow' 
lines in: (a.) Control run; (b.-d.) .AESCCHNG Mod assimilating fSCA by attempting to 
preserve SWE; (e.-f.) Enhanced Mod where new fSCA=1.0 enforced if SWE > SI; and (g.-h.) 
Runs assimilating fSCA to derived new SWE using the ADC. 

3. Components of a satellite's electro-optically sensed snow extent value (top panel) in a 
conceptual (SNOW-17)+(SAC-HT) grid square, and of below-canopy portion V that may be 
undetected (bottom panel). 

4. Example precipitation, SWE and streamflow plots (a-c): uncalibrated Carson runs for May- 
July 2005 springtime melt season. Bottom plots d-f display corresponding monthly plots of 
mean of discharge residual improvement (by assimilation) against observed streamflow. 

5 . Correlation with observations, improvement in correlation over control, percent times when 
assimilation run streamflow better than control, and corresponding % residual improvement 
(Legend for 2-letter x-axis tick label: G=Gardnerville, M=Markleeville, B=Blind, 
A=American, U=Uncalibrated, C=Calibrated). 
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Tables 


# 

Parameter 

Values/ranges [grid scaling factor] 

Units 



GU** 

GC ** 

MU or BU** 

MC** 

BC** 

AU** 

AC** 


1 

SCF 

1 

0.75-1.5 

1 

1.07-1.25 

1-1.15 

0.75-1.15 

Same* 

- 

2 

MFMAX 

0.5-1.29 [1] 

0.45-1.6 [1] 

0. 5-1.2 [1] 

0.54-1.18 [1.075] 

0.54-1.18 [1] 

0.51-1.71 [1] 

0.51-1.71 [1.5] 

mm/ °C/(6 hr) 

3 

MFMIN 

0.2-0. 5 [1] 

0.08-0.53 [1] 

0. 2-0.4 [1] 

0.08-0.53 [1] 

0.08-0.53 [1] 

0.08-0.53 [1] 

0.08-0.53 [1.5] 

mm/ °C/(6 hr) 

4 

UADJ 

0.05 

Same* 

0.05 

Same* 

Same* 

0-0.03 

Same* 

mm/mb/(6 hr) 

5 

SI 

500 

Same* 

500 

Same* 

Same* 

999 

Same* 

mm 

6 

to 

16 

ADC 

0.05,0.15,0.29, 

0.41,0.51,0.60, 

0.65,0.68,0.72, 

0.76,1.0 

Same* 

Same as GU 

Same* 

Same* 

0.0,0.1,0.2,03, 
0.4,0.5,0.6,0.7, 
0.8, 0.9, 1.0 

Same* 

- 

17 

NMF 

0.15 

Same* 

0.15 

Same* 

Same* 

0.15 

Same* 

mm/ °C/(6 hr) 

18 

TIPM 

0.06 

0.05-0.1 

0.06 

0.05-0.1 

0.05-0.1 

0. 1-0.2 

Same* 

1/(6 hr) 

19 

PXTEMP 

2 

Same* 

2 

Same* 

Same* 

2 

Same* 

°C 

20 

MBASE 

0 

Same* 

0 

Same* 

Same* 

0 

Same* 

°c 

21 

PLWHC 

0.03 

Same* 

0.03 

Same* 

Same* 

0.05-0.3 

Same* 

- 

22 

DAYGM 

0.2 

0. 1-0.3 

0.2 

0. 1-0.3 

0. 1-0.3 

0.3 

Same* 

mm/day 


* Same as respective uncalibrated value/s 

** G=Gardnerville, M=Markleeville, B=Blind, A=American, U=Uncalibrated, C=Calibrated 


Table 1: SNOW-17 parameter values, ranges, and grid scaling factors. Parameter descriptions given in Section VI of Anderson [2006]. 
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# 

Parameter 

Values/ranges [grid scaling factor] 

Units 



GU** 

GC ** 

MU or BU** 

MC** 

BC** 

AU** 

AC** 


1 

UZTWM 

23-64 [1] 

6.48-42.56 [1.5] 

23-64 [1] 

6.48-42.56 [2.225] 

6.48-42.56 [1.5] 

41-59 [1] 

41-59 [1.218] 

mm 

2 

UZFWM 

21-54 [1] 

30.97-73.75 [2.5] 

38-54 [1] 

13.48-42.97 [5] 

44.99-63.93 [2.5] 

29-53 [1] 

29-53 [1.2] 

mm 

3 

UZK 

0.34-0.75 [1] 

0.23-0.47 [0.475] 

0.47-0.75 [1] 

0.25-0.41 [0.475] 

0.25-0.41 [0.475] 

0.38-0.51 [1] 

0.38-0.51 [0.7] 

l/day 

4 

ZPERC 

45.2-79.5 [1] 

52.2-118.4 [0.875] 

Same as GU 

52.25-118.37 [0.5] 

52.23-118.37 [0.875] 

44.2-108 [1] 

Same* 

- 

5 

REXP 

1.01-2.89 [1] 

0.99-2.55 [1] 

1.01-2.07 [1] 

0.99-2.55 [1] 

0.99-2.55 [1] 

1.6-2.49 [1] 

1.6-2.49 [0.95] 

- 

6 

LZTWM 

59-252 [1] 

28.7-122.5 [0.825] 

59-252 [1] 

24.8-106.1 [1.0875] 

28.67-122.47 [0.825] 

92-248 [1] 

92-248 [4.7] 

mm 

7 

LZFSM 

2.5-20.8 [1] 

27.69-230.42 [1] 

2.5-20.8 [1] 

27.69-230.42 [1] 

27.69-230.42 [1] 

7.9-34.9 [1] 

7.9-34.9 [8.5] 

mm 

8 

LZFPM 

15-194 [1] 

23.44-187.55 [1] 

89-194 [1] 

85.52-186.43 [1] 

85.52-186.43 [1] 

46-191 [1] 

46-191 [1.11] 

mm 

9 

LZSK 

0.14-0.26 [1] 

0.04-0.07 [1] 

0.18-0.26 [1] 

0.04-0.07 [1.39] 

0.04-0.07 [1] 

0.14-0.18 [1] 

0.14-0.18 [0.5] 

l/day 

10 

LZPK 

0-0.03 [1] 

0-0.01 [1] 

0-0.03 [1] 

0-0.005 [1] 

0-0.005 [1] 

0-0.03 [1] 

0-0.03 [0.29] 

l/day 

11 

PFREE 

0.01-0.46 [1] 

0.02-0.45 [2] 

0.01-0.18 [1] 

0.03-0.45 [1.89] 

0.03-0.45 [2] 

0.08-0.3 [1] 

0.08-0.3 [0.51] 

- 

12 

PCTIM 

0 

Same* 

0 

Same* 

Same* 

0.003 

Same* 

- 

13 

ADIMP 

0 

0-0.2 

0 

0-0.2 

0-0.2 

0-0.25 

Same* 

- 

14 

RIVA 

0 

Same* 

0 

Same* 

Same* 

0.001 

Same* 

- 

15 

SIDE 

0 

Same* 

0 

Same* 

Same* 

0 

Same* 

- 

16 

RSERV 

0.3 

Same* 

0.3 

Same* 

Same* 

0.3 

Same* 

- 


* Same as respective uncalibrated value/s 

** G=Gardnerville, M=Markleeville, B=Blind, A=American, U=Uncalibrated, C=Calibrated 


Table 2: SAC-HT parameter values, ranges, and grid scaling factors. Parameter descriptions given in Table I of Koren et al. [2004]. 
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# 

Parameter 

Values/ranges [grid scaling factor] 

Units 



GU** 

GC ** 

MU or BU** 

MC** 

BC** 

AU** 

AC** 


1 

SLOPH 

0.13-0.5 

Same* 

0.22-0.5 

Same* 

Same* 

0.09-0.3 

Same* 

- 

2 

DS 

2.5 

Same* 

2.5 

Same* 

Same* 

2.5 

Same* 

1/km 

3 

ROUGH 

0.15 

Same* 

0.15 

Same* 

Same* 

0.15 

Same* 

- 

4 

QOCHN 

0.1-1.59 

Same* 

0.11-1.34 

Same* 

Same* 

0.72-2.86 [1] 

0.7-2.86 [1.75] 

m**3/s/m**(2*QMCHN) 

5 

QMCHN 

1.79 

Same* 

1.79 

Same* 

Same* 

1.47 [1] 

1.47 [0.8] 

- 


* Same as respective uncalibrated value/s 

** G=Gardnerville, M=Markleeville, B=Blind, A=American, U=Uncalibrated, C=Calibrated 


Table 3: Routing parameter values, ranges, and grid scaling factors. Parameter descriptions given in Table 3.2 of the HL-RDHM user 

manual at http://www.cbrfc.noaa.gov/present/rdhm/RDHM_3_0_0_User_Manual.pdf 
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Run 

Run description 

Short description 

A 

Operational .AESCCHNG Mod assimilating 
MODIS snow extents 

MODIS: Mod 

B 

Enhanced Mod assimilating MODIS without 
allowance for snow under forest canopy 

MODIS: No 
UnsensedSnow 

C 

Enhanced Mod assimilating MODIS with 
allowance for snow under forest canopy 

MODIS: 

UnsensedSnow 

D 

MODIS assimilation updating SWE with allowance 
for snow under forest canopy 

MODIS To SWE: 
UnsensedSnow 

Z 

Calibrated control run (for calculating improvement 
against uncalibrated control) 

- 


Table 4: Descriptions and notations used for assimilation runs (and calibrated 
control run compared against uncalibrated control) 
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Runs: => 

A 

B 

C 

D 

z 

A 

B 

C 

D 


Outlet 

Uncalibrated 


Calibrated 

RIw (%) 

Gardnerville 

26 

-18 

-5 

42 

67 

-1113 

-57 

-16 

10 

Markleeville 

37 

-18 

-7 

48 

67 

-1288 

-5 

-7 

0 

Blind 

37 

-18 

-7 

48 

66 

-1046 

-75 

-3 

10 

American 

-403 

-3 

-5 

25 

-12 

26 

-1 

-7 

34 

Iw (%) 

Gardnerville 

5 

-4 

-1 

5 

13 

-75 

-4 

-1 

1 

Markleeville 

8 

-4 

-1 

70 

14 

-88 

0 

0 

0 

Blind 

8 

-4 

-1 

10 

14 

-73 

-5 

0 

1 

American 

-153 

-1 

-2 

10 

-4 

11 

0 

-3 

15 

RI C (%) 

Gardnerville 

78 

-37 

-8 

29 

25 

45 

-63 

-1 

12 

Markleeville 

73 

-40 

-9 

27 

20 

48 

-55 

0 

7 

Blind 

73 

-40 

-9 

27 

24 

46 

-63 

-2 

11 

American 

-3 

-242 

-21 

-3 

-234 

64 

-42 

-21 

49 

Ic (days) 

Gardnerville 

10.1 

-4.8 

-1.0 

3.7 

3.2 

4.4 

-6.1 

-0.1 

1.1 

Markleeville 

9.1 

-5.0 

-1.1 

3.4 

2.5 

4.8 

-5.5 

0.0 

0.7 

Blind 

9.1 

-5.0 

-1.1 

3.4 

3.0 

4.4 

-6.0 

-0.2 

1.1 

American 

-0.1 

-8.6 

-0.7 

-0.1 

-8.3 

7.6 

-4.9 

-2.5 

5.8 

RIy (%) 
streamflow 

Gardnerville 

-122 

-159 

7 

72 

71 

-462 

-558 

23 

-100 

Markleeville 

-86 

-111 

8 

64 

90 

-1277 

-1458 

4 

-249 

Blind 

-86 

-111 

8 

64 

83 

-773 

-881 

-14 

-64 

American 

53 

81 

5 

1 

81 

-72 

-325 

22 

-45 

Iy (%) 

streamflow 

Gardnerville 

-26 

-34 

2 

76 

15 

-29 

-35 

1 

-6 

Markleeville 

-25 

-32 

2 

18 

26 

-36 

-41 

0 

-7 

Blind 

-25 

-32 

2 

18 

24 

-38 

-43 

-1 

-3 

American 

14 

22 

1 

0 

22 

-4 

-17 

1 

-2 

RIy (%) 
flood 

Gardnerville 

-99 

-92 

-1 

19 

28 

-185 

-155 

-2 

7 

Markleeville 

-150 

-132 

3 

33 

29 

-255 

-161 

-1 

18 

Blind 

-150 

-132 

3 

33 

16 

-238 

-183 

-4 

10 

American 

0 

8 

3 

0 

81 

-2 

16 

3 

0 

Iy (%) 

flood 

Gardnerville 

-27 

-25 

0 

5 

7 

-36 

-30 

0 

1 

Markleeville 

-27 

-24 

0 

6 

5 

-33 

-21 

0 

2 

Blind 

-27 

-24 

0 

6 

3 

-36 

-28 

-1 

2 

American 

0 

6 

2 

0 

60 

0 

2 

0 

0 

Rip (%) 

Gardnerville 

-110 

-104 

-2 

26 

12 

-159 

-132 

-2 

5 

Markleeville 

-126 

-111 

2 

34 

32 

-216 

-131 

-2 

10 

Blind 

-126 

-111 

2 

34 

19 

-201 

-152 

-3 

7 

American 

2 

6 

5 

0 

78 

6 

14 

3 

2 

I P (%) 

Gardnerville 

-29 

-28 

-1 

7 

3 

-37 

-31 

0 

1 

Markleeville 

-28 

-25 

0 

8 

7 

-33 

-20 

0 

2 

Blind 

-28 

-25 

0 

8 

4 

-36 

-28 

-1 

1 

American 

2 

8 

6 

0 

95 

2 

4 

1 

1 


Table 5: Performance improvement values (see Table 4 for run notations). Bold 
values are (assumed) significant improvements, italicized values are individually 


mentioned in section 6 and 7. 





Mean of Disch. Resid. o, v 

Impr. over Control (rrP/s) Volumetric discharge (m /s) Volume (mm) Volume (mm) 


(a) Hourly precipitation (Rain shown +ve, snow -ve) 



May 1 


100 


50 
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40 


20 


-20 


-40 


(d) May 


(b) Total Water Equivalent 



SNODAS Mean 
Control 

[Z] Calibrated Control 

[B] MODIS: No UnsensedSnow 

[C] MODIS: UnsensedSnow 

[D] MODIS To SWE: UnsensedSnow 


Jun 1 


Jul 1 


(c) Streamflow 



Obs 
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1 ^ n . . 

ll 

.Jo. 
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